邏輯斯回歸的限制(1)

下面是一個能夠對任意2維分布做邏輯斯迴歸運算並分類的函數，請大家試試看不同分布的狀況下邏輯斯迴歸的分類能力。

LOGISTIC_PLOT <- function (x1, x2, y, fomula) {
  
  require(scales)
  require(plot3D)
  
  model = glm(fomula, family = 'binomial')
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  z_matrix = sapply(x2_seq, function(x) {1/(1+exp(-predict(model, data.frame(x1 = x1_seq, x2 = x))))})
  
  
  image2D(z = z_matrix,
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  
}

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = 1 + 0.5 * x1 + 2 * x2
y = lr1 > 0

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2)

邏輯斯回歸的限制(2)

讓我們試試看不同分布的效果吧：

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = 1 + 0.5 * x1 + 0.7 * x2 + 3 * x1 * x2
y = lr1 > 0

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2)

這個情境下，似乎並不是分割得非常好，我們可以透過在迴歸式中加入交互作用來解決：

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2 + x1:x2)

邏輯斯回歸的限制(3)

再換個分布吧，我們同樣嘗試具有交互作用項的邏輯斯回歸：

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = - 0.5 + 0.5 * x1^2 + 0.3 * x2^2 + 0.4 * x1 * x2
y = lr1 > 0

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2 + x1:x2)

看來加入交互作用項並非萬解，其實這個情境下還要加入平方項：

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2 + x1:x2 + poly(x1, 2) + poly(x2, 2))

神經網路介紹(1)

線性迴歸及邏輯斯回歸最大的問題在於其是線性預測方法，他的預測式限制了他的發展潛力。當遇到「線性不可分割」的問題時，他們將需要「特徵工程」的協助才能做出完美預測。

– 線性迴歸

\[\hat{y} = f(x) = b_{0} + b_{1}x_1 + b_{2}x_2\]

– 邏輯斯回歸

\[log(\frac{{p}}{1-p}) = b_{0} + b_{1}x_1 + b_{2}x_2\]

儘管我們剛剛透過「特徵工程」做出了完美預測，但在低維空間中我們還能透過肉眼敏銳的觀察，但到高維空間中這將顯得不現實，我們迫切的需要一種方法能夠幫助我們「自動」做特徵工程。

神經網路介紹(2)

這個問題從1950年代後吸引眾多學者研究，而解決的方法其實很多，而效果各異。為了更好的與後續深度學習的課程做結合，我們這裡從David Rumelhart、Geoffrey Hinton、Ronald Williams在1986年的研究開始講起（事實上這並非當時的最佳解法）。

F2_1

從電腦剛開始發展的時候，我們就已經開始思考為什麼辨認貓狗對於人類來說這麼容易的事情，讓電腦做起來卻如此困難?

F2_2

這時候生物學家開始踏入了這個領域，如果我們能了解大腦是如何運作的，那是否能模仿其結構進行分類?

神經網路介紹(3)

大腦是神經系統的一部分，而他是由多個神經元互相結合而成的。(下圖為Ramón y Cajal在1905年所畫下的神經細胞樣貌)

F2_3

– 神經細胞的構造如下，不論是何種神經元皆可分成：接收區、觸發區、傳導區和輸出區。

F2_4

神經細胞傳導訊息的過程是透過電位的變化，在接收區會接收到一個電位改變的訊號，再交給突觸來處理。

F2_5

透過樹突(dendrite)能接收上一個神經元的訊息，而有些會在接收訊息後產生抑制性作用，有些會產生興奮性作用，然後這些訊號再透過神經元整合，之後再透過軸突(axon)將訊號傳導出去。
我們根據這樣的生物學知識，開始來用電腦模擬一個簡單的神經元。

神經網路介紹(4)

第一代的人工智慧是1958年由Frank Rosenblatt所發展的感知機(perceptron)，此時還稱不上是神經「網路」，他的結構如下：

F2_6

讓我們用數學語言來描述他：

\[ \begin{align} \mbox{weighted sum} & = w_{0} + w_{1}x_1 + w_{2}x_2 + \dots \\ \hat{y} & = step(\mbox{weighted sum}) \end{align} \]

仔細一看，他跟邏輯斯回歸實在是很像，差別只是在output並非是一個機率，而是直接決定是TRUE還是FALSE

– 既然是邏輯斯回歸的簡單版，那就不用談了，我們已經很清楚邏輯斯回歸的極限了。

神經網路介紹(5)

顯然大腦的構成不是僅僅使用一個神經元，既然邏輯斯回歸的構造與神經元很像，所以我們試著把多個邏輯斯回歸結合在一起吧！

– 我們不要拿太複雜的結構，就用下面這個最簡單的結構，這被稱為多層感知機(Multilayer Perceptron, MLP)

F2_7

讓我們用數學語言來描述他，為了簡化表達此一預測方程，我們將大量的使用函數組來描述：

線性預測函數L：

\[L^k(x_1, x_2) = w_{0}^k + w_{1}^kx_1 + w_{2}^kx_2\]

邏輯斯轉換函數S：

\[ \begin{align} S(x) & = \frac{{1}}{1+e^{-x}} \end{align} \]

多層感知機預測函數：

\[ \begin{align} h_1 & = S(L^1(x_1, x_2)) \\ h_2 & = S(L^2(x_1, x_2)) \\ o & = S(L^3(h_1, h_2)) \end{align} \]

反向傳播演算法(1)

有了多層感知機的預測函數之後，我們需要再給他一個損失函數再求解。而剛剛那篇1986年的Nature文章就是結合連鎖率與相關微分技術用於多層感知機，並將他命名為反向傳播演算法(backpropagation algorithm)。在這裡我們使用的是交叉熵函數作為損失函數，我們把求解目標式寫下來：

\[ \begin{align} h_1 & = S(L^1(x_1, x_2)) \\ h_2 & = S(L^2(x_1, x_2)) \\ o & = S(L^3(h_1, h_2)) \\ loss & = CE(y, o) = \frac{{1}}{n}\sum \limits_{i=1}^{n} -\left(y_{i} \cdot log(o_{i}) + (1-y_{i}) \cdot log(1-o_{i})\right) \end{align} \]

我們的目標很明確，找出一組\(w_{0}^1\)、\(w_{1}^1\)、\(w_{2}^1\)、\(w_{0}^2\)、\(w_{1}^2\)、\(w_{2}^2\)、\(w_{0}^3\)、\(w_{1}^3\)、\(w_{2}^3\)，使的整體的loss最小化。我們一樣是使用梯度下降法進行求解，而其實就是要求得上述所有係數的偏導函數。

– 其實這種複雜的微分是有規律可循的，我們把幾個重要的微分工具寫出來：

連鎖率

\[\frac{\partial}{\partial x}h(x) = \frac{\partial}{\partial x}f(g(x)) = \frac{\partial}{\partial g(x)}f(g(x)) \cdot\frac{\partial}{\partial x}g(x)\]

邏輯斯轉換函數S微分(第一課有詳細解法)

\[ \begin{align} \frac{\partial}{\partial x}S(x) & = S(x)(1-S(x)) \end{align} \]

單一交叉熵函數的微分(第一課有詳細解法)

\[ \begin{align} \frac{\partial}{\partial p_i}CE(y_i, p_i) & = \frac{p_i - y_i}{p_i(1-p_i)} \end{align} \]

反向傳播演算法(2)

下列是\(w_{0}^1\)、\(w_{1}^1\)、\(w_{2}^1\)、\(w_{0}^2\)、\(w_{1}^2\)、\(w_{2}^2\)、\(w_{0}^3\)、\(w_{1}^3\)、\(w_{2}^3\)各自的偏微分：

– 先求離結果比較近的部分(比較深層)

\(w_{0}^3\)的偏導函數：

\[ \begin{align} \frac{\partial}{\partial w_{0}^3}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{\partial}{\partial w_{0}^3}CE(y_i, o_i)\right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{\partial}{\partial o_i}CE(y_i, o_i) \cdot \frac{\partial}{\partial w_{0}^3}S(L^3(h_{1i}, h_{2i})) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{o_i-y_i}{o_i(1-o_i)} \cdot \frac{\partial}{\partial L^3(h_{1i}, h_{2i})}S(L^3(h_{1i}, h_{2i})) \cdot \frac{\partial}{\partial w_{0}^3}L^3(h_{1i}, h_{2i}) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{o_i-y_i}{o_i(1-o_i)} \cdot S(L^3(h_{1i}, h_{2i})) \cdot (1 - S(L^3(h_{1i}, h_{2i}))) \cdot \frac{\partial}{\partial w_{0}^3} \left( w_{0}^3 + w_{1}^3h_{1i} + w_{2}^3h_{2i} \right) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{o_i-y_i}{o_i(1-o_i)} \cdot o_i(1-o_i) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \end{align} \]

\(w_{1}^3\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{1}^3}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot h_{1i} \end{align} \]

\(w_{2}^3\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{2}^3}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n}\left( o_i-y_i \right) \cdot h_{2i} \end{align} \]

– 再求離結果比較遠的部分(比較淺層)

\(w_{0}^2\)的偏導函數：

\[ \begin{align} \frac{\partial}{\partial w_{0}^2}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot \frac{\partial}{\partial w_{0}^2} \left( w_{0}^3 + w_{1}^3h_{1i} + w_{2}^3h_{2i} \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3\cdot \frac{\partial}{\partial w_{0}^2} \left( h_{2i} \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3\cdot \frac{\partial}{\partial w_{0}^2} \left( S(L^2(x_{1i}, x_{2i})) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot \frac{\partial}{\partial L^2(x_{1i}, x_{2i})} S(L^2(x_{1i}, x_{2i})) \cdot \frac{\partial}{\partial w_{0}^2} L^2(x_{1i}, x_{2i}) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot h_{2i} (1 - h_{2i}) \cdot \frac{\partial}{\partial w_{0}^2} (w_{0}^2 + w_{1}^2x_{1i} + w_{2}^2x_{2i}) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot h_{2i} (1 - h_{2i}) \end{align} \]

\(w_{1}^2\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{1}^2}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot h_{2i} (1 - h_{2i}) \cdot x_{1i} \end{align} \]

\(w_{2}^2\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{2}^2}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot h_{2i} (1 - h_{2i}) \cdot x_{2i} \end{align} \]

\(w_{0}^1\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{0}^1}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{1}^3 \cdot h_{1i} (1 - h_{1i}) \end{align} \]

\(w_{1}^1\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{1}^1}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{1}^3 \cdot h_{1i} (1 - h_{1i}) \cdot x_{1i} \end{align} \]

\(w_{2}^1\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{2}^1}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{1}^3 \cdot h_{1i} (1 - h_{1i}) \cdot x_{2i} \end{align} \]

練習1：利用梯度下降法獲得一個多層感知機

怕你看不懂公式，幫妳寫成公式的樣子

#Sample generation

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = - 1.5 + x1^2 + x2^2
y = lr1 > 0 + 0L

#Forward

S.fun = function (x, eps = 1e-5) {
  S = 1/(1 + exp(-x))
  S[S < eps] = eps
  S[S > 1 - eps] = 1 - eps
  return(S)
}

h1.fun = function (w10, w11, w12, x1 = x1, x2 = x2) {
  L1 = w10 + w11 * x1 + w12 * x2
  return(S.fun(L1))
}

h2.fun = function (w20, w21, w22, x1 = x1, x2 = x2) {
  L2 = w20 + w21 * x1 + w22 * x2
  return(S.fun(L2))
}

o.fun = function (w30, w31, w32, h1, h2) {
  L3 = w30 + w31 * h1 + w32 * h2
  return(S.fun(L3))
}

loss.fun = function (o, y = y) {
  loss = -1/length(y) * sum(y * log(o) + (1 - y) * log(1 - o))
  return(loss)
}

#Backward

differential.fun.w30 = function(o, y = y) {
  return(-1/length(y)*sum(y-o))
}

differential.fun.w31 = function(o, h1, y = y) {
  return(-1/length(y)*sum((y-o)*h1))
}

differential.fun.w32 = function(o, h2, y = y) {
  return(-1/length(y)*sum((y-o)*h2))
}

differential.fun.w20 = function(o, h2, w32, y = y) {
  return(-1/length(y)*sum((y-o)*w32*h2*(1-h2)))
}

differential.fun.w21 = function(o, h2, w32, y = y, x1 = x1) {
  return(-1/length(y)*sum((y-o)*w32*h2*(1-h2)*x1))
}

differential.fun.w22 = function(o, h2, w32, y = y, x2 = x2) {
  return(-1/length(y)*sum((y-o)*w32*h2*(1-h2)*x2))
}

differential.fun.w10 = function(o, h1, w31, y = y) {
  return(-1/length(y)*sum((y-o)*w31*h1*(1-h1)))
}

differential.fun.w11 = function(o, h1, w31, y = y, x1 = x1) {
  return(-1/length(y)*sum((y-o)*w31*h1*(1-h1)*x1))
}

differential.fun.w12 = function(o, h1, w31, y = y, x2 = x2) {
  return(-1/length(y)*sum((y-o)*w31*h1*(1-h1)*x2))
}

練習1答案

num.iteration = 10000
lr = 0.1
W_matrix = matrix(0, nrow = num.iteration + 1, ncol = 9)
loss_seq = rep(0, num.iteration)
colnames(W_matrix) = c('w10', 'w11', 'w12', 'w20', 'w21', 'w22', 'w30', 'w31', 'w32')

#Start random values
W_matrix[1,] = rnorm(9, sd = 1) 

for (i in 2:(num.iteration+1)) {
  
  #Forward
  
  current_H1 = h1.fun(w10 = W_matrix[i-1,1], w11 = W_matrix[i-1,2], w12 = W_matrix[i-1,3],
                      x1 = x1, x2 = x2)
  
  current_H2 = h2.fun(w20 = W_matrix[i-1,4], w21 = W_matrix[i-1,5], w22 = W_matrix[i-1,6],
                      x1 = x1, x2 = x2)
  
  current_O = o.fun(w30 = W_matrix[i-1,7], w31 = W_matrix[i-1,8], w32 = W_matrix[i-1,9],
                    h1 = current_H1, h2 = current_H2)
  
  loss_seq[i-1] = loss.fun(o = current_O, y = y)
  
  #Backward
  
  W_matrix[i,1] = W_matrix[i-1,1] - lr * differential.fun.w10(o = current_O, h1 = current_H1,
                                       w31 = W_matrix[i-1,8], y = y)
  
  W_matrix[i,2] = W_matrix[i-1,2] - lr * differential.fun.w11(o = current_O, h1 = current_H1,
                                       w31 = W_matrix[i-1,8], y = y, x1 = x1)
  
  W_matrix[i,3] = W_matrix[i-1,3] - lr * differential.fun.w12(o = current_O, h1 = current_H1,
                                       w31 = W_matrix[i-1,8], y = y, x2 = x2)
  
  W_matrix[i,4] = W_matrix[i-1,4] - lr * differential.fun.w20(o = current_O, h2 = current_H2,
                                       w32 = W_matrix[i-1,9], y = y)
    
  W_matrix[i,5] = W_matrix[i-1,5] - lr * differential.fun.w21(o = current_O, h2 = current_H2,
                                       w32 = W_matrix[i-1,9], y = y, x1 = x1)
  
  W_matrix[i,6] = W_matrix[i-1,6] - lr * differential.fun.w22(o = current_O, h2 = current_H2,
                                       w32 = W_matrix[i-1,9], y = y, x2 = x2)
    
  W_matrix[i,7] = W_matrix[i-1,7] - lr * differential.fun.w30(o = current_O, y = y)
    
  W_matrix[i,8] = W_matrix[i-1,8] - lr * differential.fun.w31(o = current_O, h1 = current_H1, y = y)
  
  W_matrix[i,9] = W_matrix[i-1,9] - lr * differential.fun.w32(o = current_O, h2 = current_H2, y = y)
  
}

圖像化

require(scales)
require(plot3D)
  
x1_seq = seq(min(x1), max(x1), length.out = 100)
x2_seq = seq(min(x2), max(x2), length.out = 100)

pre_func = function (x1, x2, W_list = W_matrix[nrow(W_matrix),]) {
  H1 = h1.fun(w10 = W_list[1], w11 = W_list[2], w12 = W_list[3], x1 = x1, x2 = x2)
  H2 = h2.fun(w20 = W_list[4], w21 = W_list[5], w22 = W_list[6], x1 = x1, x2 = x2)
  O = o.fun(w30 = W_list[7], w31 = W_list[8], w32 = W_list[9], h1 = H1, h2 = H2)
  return(O)
}

z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  
image2D(z = z_matrix,
        x = x1_seq, xlab = 'x1',
        y = x2_seq, ylab = 'x2',
        shade = 0.2, rasterImage = TRUE,
        col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))

points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)

這真是太有趣了，組合3個邏輯斯回歸讓分割邊界產生曲線，讓我們觀察一下\(loss\)下降的軌跡

plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')

矩陣化多層感知機(1)

透過多個人工神經元(邏輯斯迴歸)達到的效果令人驚嘆，但每次增加神經元時都將導致微分過程的變化，從而導致其難以擴展。

– 透過矩陣能夠將多層感知機變成較為簡易的公式組合(以下為個人層級的方程式)：

線性預測函數L之矩陣式(d為輸出維度)：

\[ \begin{align} L^k_d(X, W^k_d) & = XW^k_d \\ X & = \begin{pmatrix} x_{1,1} & x_{1,2} & \cdots & x_{1,m} \end{pmatrix} \\ W^k_d & = \begin{pmatrix} w^k_{0,1} & w^k_{0,2} & \cdots & w^k_{0,d} \\ w^k_{1,1} & w^k_{1,2} & \cdots & w^k_{1,d} \\ w^k_{2,1} & w^k_{2,2} & \cdots & w^k_{2,d} \\ \vdots & \vdots & \ddots & \vdots \\ w^k_{m,1} & w^k_{m,2} & \cdots & w^k_{m,d} \end{pmatrix} \\ \frac{\partial}{\partial W^k_d}L^k_d(X) & = \begin{pmatrix} X^T & & X^T & \cdots & X^T \\ \end{pmatrix} \mbox{ [repeat } d \mbox{ times]} \end{align} \]

邏輯斯轉換函數S：

\[ \begin{align} S(x) & = \frac{{1}}{1+e^{-x}} \\ \frac{\partial}{\partial x}S(x) & = S(x)(1-S(x)) \end{align} \]

多層感知機預測函數之矩陣式，矩陣上標E代表該矩陣的第一欄已被填滿1(這是剛剛的單一隱藏層網路，其隱藏層中神經元數目為d，在剛剛的例子中d=2)：

\[ \begin{align} l_1 & = L^1_d(x^E,W^1_d) \\ h_1 & = S(l_1) \\ l_2 & = L^2_1(h_1^E,W^2_1) \\ o & = S(l_2) \end{align} \]

矩陣化多層感知機(2)

我們同樣的需要解矩陣的微分方程，求解目標式如下：

\[ \begin{align} l_1 & = L^1_d(x^E,W^1_d) \\ h_1 & = S(l_1) \\ l_2 & = L^2_1(h_1^E,W^2_1) \\ o & = S(l_2) \\ loss & = CE(y, o) = -\left(y \cdot log(o) + (1-y) \cdot log(1-o)\right) \end{align} \]

讓我們分別對這裡所有元素進行微分(符號\(\otimes\)為按元素相乘，符號\(\bullet\)為按矩陣乘法)：

\[ \begin{align} grad.o & = \frac{\partial}{\partial o}loss = \frac{o-y}{o(1-o)} \\ grad.l_2 & = \frac{\partial}{\partial l_2}loss = grad.o \otimes \frac{\partial}{\partial l_2}o= o-y \\ grad.W^2_1 & = \frac{\partial}{\partial W^2_1}loss = grad.l_2 \otimes \frac{\partial}{\partial W^2_1}l_2 = \frac{1}{n} \otimes (h_1^E)^T \bullet grad.l_2\\ grad.h_1^E & = \frac{\partial}{\partial h_1^E}loss = grad.l_2 \otimes \frac{\partial}{\partial h_1^E}l_2 = grad.l_2 \bullet (W^2_1)^T \\ grad.l_1 & = \frac{\partial}{\partial l_1}loss = grad.h_1 \otimes \frac{\partial}{\partial l_1}h_1 = grad.h_1 \otimes h_1 \otimes (1-h_1) \\ grad.W^1_d & = \frac{\partial}{\partial W^1_d}loss = grad.l_1 \otimes \frac{\partial}{\partial W^1_d}l_1 = \frac{1}{n} \otimes (x^E)^T \bullet grad.l_1 \end{align} \]

練習2：利用矩陣式求解並試著增加隱藏層的數目

試著實現矩陣化的MLP吧！

– 我們用程式碼描述上述式子：

#Sample generation

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = - 1.5 + x1^2 + x2^2
y = lr1 > 0 + 0L

#Forward

S.fun = function (x, eps = 1e-5) {
  S = 1/(1 + exp(-x))
  S[S < eps] = eps
  S[S > 1 - eps] = 1 - eps
  return(S)
}

L.fun = function (X, W) {
  X.E = cbind(1, X)
  L = X.E %*% W
  return(L)
}

CE.fun = function (o, y, eps = 1e-5) {
  loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
  return(loss)
}

#Backward

grad_o.fun = function (o, y) {
  return((o - y)/(o*(1-o)))
}

grad_l2.fun = function (grad_o, o) {
  return(grad_o*(o*(1-o)))
}

grad_W2.fun = function (grad_l2, h1) {
  h1.E = cbind(1, h1)
  return(t(h1.E) %*% grad_l2/nrow(h1))
}

grad_h1.fun = function (grad_l2, W2) {
  return(grad_l2 %*% t(W2[-1,]))
}

grad_l1.fun = function (grad_h1, h1) {
  return(grad_h1*(h1*(1-h1)))
}

grad_W1.fun = function (grad_l1, x) {
  x.E = cbind(1, x)
  return(t(x.E) %*% grad_l1/nrow(x))
}

練習2答案

這是一個MLP訓練函數，由於已經矩陣化了，我們可以很方便的指定隱藏層的神經元數量(num.hidden)

MLP_Trainer = function (num.iteration = 500, num.hidden = 2, lr = 0.1, x1 = x1, x2 = x2, y = y) {
  
  #Functions
  
  #Forward
  
  S.fun = function (x, eps = 1e-5) {
    S = 1/(1 + exp(-x))
    S[S < eps] = eps
    S[S > 1 - eps] = 1 - eps
    return(S)
  }
  
  L.fun = function (X, W) {
    X.E = cbind(1, X)
    L = X.E %*% W
    return(L)
  }
  
  CE.fun = function (o, y, eps = 1e-5) {
    loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
    return(loss)
  }
  
  #Backward
  
  grad_o.fun = function (o, y) {
    return((o - y)/(o*(1-o)))
  }
  
  grad_l2.fun = function (grad_o, o) {
    return(grad_o*(o*(1-o)))
  }
  
  grad_W2.fun = function (grad_l2, h1) {
    h1.E = cbind(1, h1)
    return(t(h1.E) %*% grad_l2/nrow(h1))
  }
  
  grad_h1.fun = function (grad_l2, W2) {
    return(grad_l2 %*% t(W2[-1,]))
  }
  
  grad_l1.fun = function (grad_h1, h1) {
    return(grad_h1*(h1*(1-h1)))
  }
  
  grad_W1.fun = function (grad_l1, x) {
    x.E = cbind(1, x)
    return(t(x.E) %*% grad_l1/nrow(x))
  }
  
  #Caculating
  
  X_matrix = cbind(x1, x2)
  
  W1_list = list()
  W2_list = list()
  loss_seq = rep(0, num.iteration)
  
  #Start random values
  
  W1_list[[1]] = matrix(rnorm(3*num.hidden, sd = 1), nrow = 3, ncol = num.hidden)
  W2_list[[1]] = matrix(rnorm(num.hidden + 1, sd = 1), nrow = num.hidden + 1, ncol = 1)
  
  for (i in 2:(num.iteration+1)) {
    
    #Forward
    
    current_l1 = L.fun(X = X_matrix, W = W1_list[[i - 1]])
    current_h1 = S.fun(x = current_l1)
    current_l2 = L.fun(X = current_h1, W = W2_list[[i - 1]])
    current_o = S.fun(x = current_l2)
    loss_seq[i-1] = CE.fun(o = current_o, y = y, eps = 1e-5)
    
    #Backward
    
    current_grad_o = grad_o.fun(o = current_o, y = y)
    current_grad_l2 = grad_l2.fun(grad_o = current_grad_o, o = current_o)
    current_grad_W2 = grad_W2.fun(grad_l2 = current_grad_l2, h1 = current_h1)
    current_grad_h1 = grad_h1.fun(grad_l2 = current_grad_l2, W2 = W2_list[[i - 1]])
    current_grad_l1 = grad_l1.fun(grad_h1 = current_grad_h1, h1 = current_h1)
    current_grad_W1 = grad_W1.fun(grad_l1 = current_grad_l1, x = X_matrix)
    
    W2_list[[i]] = W2_list[[i-1]] - lr * current_grad_W2
    W1_list[[i]] = W1_list[[i-1]] - lr * current_grad_W1
    
  }
  
  require(scales)
  require(plot3D)
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  pre_func = function (x1, x2, W1 = W1_list[[length(W1_list)]], W2 = W2_list[[length(W2_list)]]) {
    new_X = cbind(x1, x2)
    O = S.fun(x = L.fun(X = S.fun(x = L.fun(X = new_X, W = W1)), W = W2))
    return(O)
  }
  
  z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  par(mfrow = c(1, 2))
  
  image2D(z = z_matrix,
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  
  plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')
  
}

看結果

MLP_Trainer(num.iteration = 10000, num.hidden = 2, lr = 0.1, x1 = x1, x2 = x2, y = y)

現在讓我們試試看多給他幾個神經元，是不是又切的更好了?

MLP_Trainer(num.iteration = 10000, num.hidden = 5, lr = 0.1, x1 = x1, x2 = x2, y = y)

訓練樣本與測試樣本(1)

有了MLP之後實在是太方便了，現在我們只要透過簡單的「增加神經元」(不就改個參數?)，就能做出任意形狀的分類任務。
這時候我們要思考的是，假設我們為了讓MLP能夠完美的分割出所有樣本，從而讓他學習出離譜的邊界，那會不會造成未來在應用時造成泛用性不佳。

– 舉例來說，假設我們希望做出一個良好的分類器區分下圖的藍點與紅點，你是否同意黑線相較於綠線是一條更好的區分邊界?

F2_8

訓練樣本與測試樣本(2)

現在讓我們試著為樣本加點雜訊，並且看看剛剛那個MLP分類器的效果：

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = - 1.5 + x1^2 + x2^2 + rnorm(50)
y = lr1 > 0 + 0L

MLP_Trainer(num.iteration = 10000, num.hidden = 5, lr = 0.1, x1 = x1, x2 = x2, y = y)

好像有許多樣本被錯誤分類了，讓我們試試看用100個神經元來解決事情吧！

MLP_Trainer(num.iteration = 10000, num.hidden = 100, lr = 0.1, x1 = x1, x2 = x2, y = y)

準確性的確有所提升，但透過這種奇怪的邊界所得到的提升真的有比較好嗎?

訓練樣本與測試樣本(3)

在這裡我們必須給一個觀念，在強大的「非線性分類器」(或稱作人工智慧)，只要你願意給他足夠的「記憶體」(或稱作腦容量)，你絕對可以把所有的訓練樣本都正確分類，但面對未見過的樣本就不見得有同樣良好的分布能力了。

– 讓我們用同樣的分布產生資料，並把前50個樣本作為訓練集(Training set)，後50個做測試集(Test set)

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = - 1.5 + x1^2 + x2^2 + rnorm(50)
y = lr1 > 0 + 0L

test_x1 = rnorm(50, sd = 1) 
test_x2 = rnorm(50, sd = 1) 
lr1 = - 1.5 + test_x1^2 + test_x2^2 + rnorm(50)
test_y = lr1 > 0 + 0L

我們重新寫一個MLP訓練函數，讓他能夠呈現測試集的預測結果

MLP_Trainer = function (num.iteration = 500, num.hidden = 2, lr = 0.1, x1 = x1, x2 = x2, y = y, test_x1 = NULL, test_x2 = NULL, test_y = NULL) {
  
  #Functions
  
  #Forward
  
  S.fun = function (x, eps = 1e-5) {
    S = 1/(1 + exp(-x))
    S[S < eps] = eps
    S[S > 1 - eps] = 1 - eps
    return(S)
  }
  
  L.fun = function (X, W) {
    X.E = cbind(1, X)
    L = X.E %*% W
    return(L)
  }
  
  CE.fun = function (o, y, eps = 1e-5) {
    loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
    return(loss)
  }
  
  #Backward
  
  grad_o.fun = function (o, y) {
    return((o - y)/(o*(1-o)))
  }
  
  grad_l2.fun = function (grad_o, o) {
    return(grad_o*(o*(1-o)))
  }
  
  grad_W2.fun = function (grad_l2, h1) {
    h1.E = cbind(1, h1)
    return(t(h1.E) %*% grad_l2/nrow(h1))
  }
  
  grad_h1.fun = function (grad_l2, W2) {
    return(grad_l2 %*% t(W2[-1,]))
  }
  
  grad_l1.fun = function (grad_h1, h1) {
    return(grad_h1*(h1*(1-h1)))
  }
  
  grad_W1.fun = function (grad_l1, x) {
    x.E = cbind(1, x)
    return(t(x.E) %*% grad_l1/nrow(x))
  }
  
  #Caculating
  
  X_matrix = cbind(x1, x2)
  
  W1_list = list()
  W2_list = list()
  loss_seq = rep(0, num.iteration)
  
  #Start random values
  
  W1_list[[1]] = matrix(rnorm(3*num.hidden, sd = 1), nrow = 3, ncol = num.hidden)
  W2_list[[1]] = matrix(rnorm(num.hidden + 1, sd = 1), nrow = num.hidden + 1, ncol = 1)
  
  for (i in 2:(num.iteration+1)) {
    
    #Forward
    
    current_l1 = L.fun(X = X_matrix, W = W1_list[[i - 1]])
    current_h1 = S.fun(x = current_l1)
    current_l2 = L.fun(X = current_h1, W = W2_list[[i - 1]])
    current_o = S.fun(x = current_l2)
    loss_seq[i-1] = CE.fun(o = current_o, y = y, eps = 1e-5)
    
    #Backward
    
    current_grad_o = grad_o.fun(o = current_o, y = y)
    current_grad_l2 = grad_l2.fun(grad_o = current_grad_o, o = current_o)
    current_grad_W2 = grad_W2.fun(grad_l2 = current_grad_l2, h1 = current_h1)
    current_grad_h1 = grad_h1.fun(grad_l2 = current_grad_l2, W2 = W2_list[[i - 1]])
    current_grad_l1 = grad_l1.fun(grad_h1 = current_grad_h1, h1 = current_h1)
    current_grad_W1 = grad_W1.fun(grad_l1 = current_grad_l1, x = X_matrix)
    
    W2_list[[i]] = W2_list[[i-1]] - lr * current_grad_W2
    W1_list[[i]] = W1_list[[i-1]] - lr * current_grad_W1
    
  }
  
  require(scales)
  require(plot3D)
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  pre_func = function (x1, x2, W1 = W1_list[[length(W1_list)]], W2 = W2_list[[length(W2_list)]]) {
    new_X = cbind(x1, x2)
    O = S.fun(x = L.fun(X = S.fun(x = L.fun(X = new_X, W = W1)), W = W2))
    return(O)
  }
  
  pred_y = pre_func(x1 = x1, x2 = x2)
  MAIN_TXT = paste0('Train-Acc:', formatC(mean((pred_y > 0.5) == y), 2, format = 'f'))
  if (!is.null(test_x1)) {
    pred_test_y = pre_func(x1 = test_x1, x2 = test_x2)
    MAIN_TXT = paste0(MAIN_TXT, '; Test-Acc:', formatC(mean((pred_test_y > 0.5) == test_y), 2, format = 'f'))
  }
  
  z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  par(mfrow = c(1, 2))
  
  image2D(z = z_matrix, main = MAIN_TXT,
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  if (!is.null(test_x1)) {
    points(test_x1, test_x2, col = 'black', bg = c('#C00000', '#0000C0')[(test_y + 1)], pch = 21)
  }
  
  plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')
  
}

看看結果如何：

MLP_Trainer(num.iteration = 10000, num.hidden = 5, lr = 0.1, x1 = x1, x2 = x2, y = y, test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

MLP_Trainer(num.iteration = 10000, num.hidden = 100, lr = 0.1, x1 = x1, x2 = x2, y = y, test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

練習3：如何決定最佳的神經元數量

所以我們現在知道了在訓練樣本中預測精準的模型，到了真實世界(測試樣本)的能力不見得也會一樣好。那現在的問題出來了，以剛剛我要如何決定該用多少個神經元來建構模型?

– 注意，我們需要在不使用測試集數據的狀況下獲得這樣的結果。

請你好好想一想，並且提出一個解答方案。

練習3答案

在不使用測試集的前提下，我們只能在原有的樣本中再切割出一小部分作為「驗證集(Validation set)」進行神經元數量不同時的外推性比較，以此確定神經元數目。

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = - 1.5 + x1^2 + x2^2 + rnorm(50)
y = lr1 > 0 + 0L

test_x1 = rnorm(50, sd = 1) 
test_x2 = rnorm(50, sd = 1) 
lr1 = - 1.5 + test_x1^2 + test_x2^2 + rnorm(50)
test_y = lr1 > 0 + 0L

MLP_Trainer(num.iteration = 10000, num.hidden = 5, lr = 0.1, x1 = x1[1:30], x2 = x2[1:30], y = y[1:30], test_x1 = x1[31:50], test_x2 = x2[31:50], test_y = y[31:50])

MLP_Trainer(num.iteration = 10000, num.hidden = 20, lr = 0.1, x1 = x1[1:30], x2 = x2[1:30], y = y[1:30], test_x1 = x1[31:50], test_x2 = x2[31:50], test_y = y[31:50])

結果看起來神經元數量為5是最佳的模型，讓我們真的把全訓練樣本都用神經元數量為5及20再訓練一次：

MLP_Trainer(num.iteration = 10000, num.hidden = 5, lr = 0.1, x1 = x1, x2 = x2, y = y, test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

MLP_Trainer(num.iteration = 10000, num.hidden = 20, lr = 0.1, x1 = x1, x2 = x2, y = y, test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

非線性轉換函數的問題(1)

在剛剛的眾多實驗中，敏銳的你應該發現了一個現象，也就是我們在訓練(隨著梯度找尋最低點)的過程中，損失值經常會先經過一些瓶頸再急遽下降。

– 這個問題其實我們在第一課已經遇到了，當時我們使用殘差平方和作為邏輯斯回歸的損失函數，導致偏導函數中出現了\(p(1-p)\)的部分，而我們的解決方式是透過改寫損失函數將這個部分約分掉，從而讓梯度實現平穩的下降。

讓我們看看我們目前的偏導函數，你會發現\(W^1_d\)的偏導函數中同樣存在了\(p(1-p)\)的部分，而出線的原因也很容易想到，那是因為在隱藏層與輸出層的中間我們使用了邏輯斯轉換函數進行轉換。

\[ \begin{align} grad.l_1 & = \frac{\partial}{\partial l_1}loss = grad.h_1 \otimes \frac{\partial}{\partial l_1}h_1 = grad.h_1 \otimes h_1 \otimes (1-h_1) \\ grad.W^1_d & = \frac{\partial}{\partial W^1_d}loss = grad.l_1 \otimes \frac{\partial}{\partial W^1_d}l_1 = \frac{1}{n} \otimes (x^E)^T \bullet grad.l_1 \end{align} \]

想想我們該如何解決這個問題吧！

非線性轉換函數的問題(2)

我們不太可能再次改寫損失函數達到同樣的效果，因為以剛剛的結構為例，第一層的輸出經過一次邏輯斯轉換函數的轉換，而第二層則經過了兩次。因此比較合理的做法是修改中間層的轉換函數。

– 線性整流函數(Rectified Linear Unit, ReLU)首次被提出的時間已經不可考了，但其實際開始大量被應用於中間層的轉換函數是2012年Alexnet所開始的。AlexNet可以說是現代深度神經網路的開山之作，並揭開了深度學習的熱潮。 2012年他一舉摘下了ILSVRC競賽的冠軍，並且效果大幅度超過傳統的方法，也讓錯誤率從25%降低至15%以下。

F2_9

– ReLU的函數形式：

\[ ReLU(x) = \left\{ \begin{array} -x & \mbox{ if x > 0} \\ 0 & \mbox{ otherwise} \end{array} \right. \]

我們可以想像，在一個足夠深的網路之中，中間經過了非常多次邏輯斯轉換函數的轉換，那對於最前面幾層的係數將會非常難以訓練。這是因為在訓練過程中只要每經過一次邏輯斯轉換，那就會增加了梯度消失的風險。而使用ReLU則能避免這個問題。

– 注意，這個問題叫做「梯度消失問題」(gradient vanishing problem)，這個問題是Deep Learning的經典問題，而其實問題一直就是出在這個地方。問題直到了2016年開始才有比較理想的解決方案。

非線性轉換函數的問題(3)

ReLU函數的導函數相當簡單，而他能很輕易的與我們之前的多層感知機迅速結合：

– ReLU的導函數

\[ \frac{\partial}{\partial x}ReLU(x) = \left\{ \begin{array} -1 & \mbox{ if x > 0} \\ 0 & \mbox{ otherwise} \end{array} \right. \]

– 使用ReLU作為轉換函數的MLP預測方程：

\[ \begin{align} l_1 & = L^1_d(x^E,W^1_d) \\ h_1 & = ReLU(l_1) \\ l_2 & = L^2_1(h_1^E,W^2_1) \\ o & = S(l_2) \\ loss & = CE(y, o) = -\left(y \cdot log(o) + (1-y) \cdot log(1-o)\right) \end{align} \]

– 各元素的導函數：

\[ \begin{align} grad.o & = \frac{\partial}{\partial o}loss = \frac{o-y}{o(1-o)} \\ grad.l_2 & = \frac{\partial}{\partial l_2}loss = grad.o \otimes \frac{\partial}{\partial l_2}o= o-y \\ grad.W^2_1 & = \frac{\partial}{\partial W^2_1}loss = grad.l_2 \otimes \frac{\partial}{\partial W^2_1}l_2 = \frac{{1}}{n} \otimes (h_1^E)^T \bullet grad.l_2\\ grad.h_1^E & = \frac{\partial}{\partial h_1^E}loss = grad.l_2 \otimes \frac{\partial}{\partial h_1^E}l_2 = grad.l_2 \bullet (W^2_1)^T \\ grad.l_1 & = \frac{\partial}{\partial l_1}loss = grad.h_1 \otimes \frac{\partial}{\partial l_1}h_1 = grad.h_1 \otimes \frac{\partial}{\partial l_1}ReLU(l_1) \\ grad.W^1_d & = \frac{\partial}{\partial W^1_d}loss = grad.l_1 \otimes \frac{\partial}{\partial W^1_d}l_1 = \frac{{1}}{n} \otimes (x^E)^T \bullet grad.l_1 \end{align} \]

練習4：建構使用ReLU的MLP函數

使用ReLU作為轉換函數的MLP預測函數及其偏導函數如下，請你試著建構使用ReLU的MLP求解函數，並做實驗看看其損失函數的變化情形。

#Sample generation

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
X = cbind(x1, x2)
lr1 = - 1.5 + x1^2 + x2^2
y = lr1 > 0 + 0L

#Forward

S.fun = function (x, eps = 1e-5) {
  S = 1/(1 + exp(-x))
  S[S < eps] = eps
  S[S > 1 - eps] = 1 - eps
  return(S)
}

ReLU.fun = function (x) {
  x[x < 0] <- 0
  return(x)
}

L.fun = function (X, W) {
  X.E = cbind(1, X)
  L = X.E %*% W
  return(L)
}

CE.fun = function (o, y, eps = 1e-5) {
  loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
  return(loss)
}

#Backward

grad_o.fun = function (o, y) {
  return((o - y)/(o*(1-o)))
}

grad_l2.fun = function (grad_o, o) {
  return(grad_o*(o*(1-o)))
}

grad_W2.fun = function (grad_l2, h1) {
  h1.E = cbind(1, h1)
  return(t(h1.E) %*% grad_l2/nrow(h1))
}

grad_h1.fun = function (grad_l2, W2) {
  return(grad_l2 %*% t(W2[-1,]))
}

grad_l1.fun = function (grad_h1, l1) {
  de_l1 = l1
  de_l1[de_l1<0] = 0
  de_l1[de_l1>0] = 1
  return(grad_h1*de_l1)
}

grad_W1.fun = function (grad_l1, x) {
  x.E = cbind(1, x)
  return(t(x.E) %*% grad_l1/nrow(x))
}

練習4答案

這是使用ReLU作為轉換函數的MLP預測函數及其偏導函數如下，請你試著建構使用ReLU的MLP求解函數，並做實驗看看其損失函數的變化情形。

MLP_Trainer = function (num.iteration = 500, num.hidden = 2, lr = 0.1, x1 = x1, x2 = x2, y = y, test_x1 = NULL, test_x2 = NULL, test_y = NULL) {
  
  #Functions
  
  #Forward
  
  S.fun = function (x, eps = 1e-5) {
    S = 1/(1 + exp(-x))
    S[S < eps] = eps
    S[S > 1 - eps] = 1 - eps
    return(S)
  }
  
  ReLU.fun = function (x) {
    x[x < 0] <- 0
    return(x)
  }
  
  L.fun = function (X, W) {
    X.E = cbind(1, X)
    L = X.E %*% W
    return(L)
  }
  
  CE.fun = function (o, y, eps = 1e-5) {
    loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
    return(loss)
  }
  
  #Backward
  
  grad_o.fun = function (o, y) {
    return((o - y)/(o*(1-o)))
  }
  
  grad_l2.fun = function (grad_o, o) {
    return(grad_o*(o*(1-o)))
  }
  
  grad_W2.fun = function (grad_l2, h1) {
    h1.E = cbind(1, h1)
    return(t(h1.E) %*% grad_l2/nrow(h1))
  }
  
  grad_h1.fun = function (grad_l2, W2) {
    return(grad_l2 %*% t(W2[-1,]))
  }
  
  grad_l1.fun = function (grad_h1, l1) {
    de_l1 = l1
    de_l1[de_l1<0] = 0
    de_l1[de_l1>0] = 1
    return(grad_h1*de_l1)
  }
  
  grad_W1.fun = function (grad_l1, x) {
    x.E = cbind(1, x)
    return(t(x.E) %*% grad_l1/nrow(x))
  }
  
  #Caculating
  
  X_matrix = cbind(x1, x2)
  
  W1_list = list()
  W2_list = list()
  loss_seq = rep(0, num.iteration)
  
  #Start random values
  
  W1_list[[1]] = matrix(rnorm(3*num.hidden, sd = 1), nrow = 3, ncol = num.hidden)
  W2_list[[1]] = matrix(rnorm(num.hidden + 1, sd = 1), nrow = num.hidden + 1, ncol = 1)
  
  for (i in 2:(num.iteration+1)) {
    
    #Forward
    
    current_l1 = L.fun(X = X_matrix, W = W1_list[[i - 1]])
    current_h1 = ReLU.fun(x = current_l1)
    current_l2 = L.fun(X = current_h1, W = W2_list[[i - 1]])
    current_o = S.fun(x = current_l2)
    loss_seq[i-1] = CE.fun(o = current_o, y = y, eps = 1e-5)
    
    #Backward
    
    current_grad_o = grad_o.fun(o = current_o, y = y)
    current_grad_l2 = grad_l2.fun(grad_o = current_grad_o, o = current_o)
    current_grad_W2 = grad_W2.fun(grad_l2 = current_grad_l2, h1 = current_h1)
    current_grad_h1 = grad_h1.fun(grad_l2 = current_grad_l2, W2 = W2_list[[i - 1]])
    current_grad_l1 = grad_l1.fun(grad_h1 = current_grad_h1, l1 = current_l1)
    current_grad_W1 = grad_W1.fun(grad_l1 = current_grad_l1, x = X_matrix)
    W2_list[[i]] = W2_list[[i-1]] - lr * current_grad_W2
    W1_list[[i]] = W1_list[[i-1]] - lr * current_grad_W1
    
  }
  
  require(scales)
  require(plot3D)
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  pre_func = function (x1, x2, W1 = W1_list[[length(W1_list)]], W2 = W2_list[[length(W2_list)]]) {
    new_X = cbind(x1, x2)
    O = S.fun(x = L.fun(X = ReLU.fun(x = L.fun(X = new_X, W = W1)), W = W2))
    return(O)
  }
  
  pred_y = pre_func(x1 = x1, x2 = x2)
  MAIN_TXT = paste0('Train-Acc:', formatC(mean((pred_y > 0.5) == y), 2, format = 'f'))
  if (!is.null(test_x1)) {
    pred_test_y = pre_func(x1 = test_x1, x2 = test_x2)
    MAIN_TXT = paste0(MAIN_TXT, '; Test-Acc:', formatC(mean((pred_test_y > 0.5) == test_y), 2, format = 'f'))
  }
  
  z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  par(mfrow = c(1, 2))
  
  image2D(z = z_matrix, main = MAIN_TXT,
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  if (!is.null(test_x1)) {
    points(test_x1, test_x2, col = 'black', bg = c('#C00000', '#0000C0')[(test_y + 1)], pch = 21)
  }
  
  plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')
  
}

看看不同的轉換函數對分類邊界產生的效果吧

MLP_Trainer(num.iteration = 10000, num.hidden = 5, lr = 0.1, x1 = x1, x2 = x2, y = y)

除了分類邊界的影響，你是否發現運算速度快了非常多?

結語

本次課程介紹了經典的反向傳播演算法應用於多層感知機求解的過程，並實際用其解決一般邏輯斯回歸較難解決的問題上。

– 另外與傳統的統計推論最大的不同點在於，我們目前掌握了強大的非線性分類器，而只要你願意他就能將整個既有樣本中的所有資訊都透過某種形式記憶住，從而達到極高的準確性，但其對於真實世界中新樣本的外推性就會被質疑。

最後我們留一個問題給大家帶回家思考：邏輯斯迴歸與神經元構造的相似性讓我們選擇了透過堆疊邏輯斯迴歸打造人工神經網路，但其導函數的性質讓現在深度神經網路大多使用ReLU作為各層之間的轉換函數。那為什麼我們不考慮在各層之間「不轉換」呢?

– 另外，我們是否有可能使用Step函數(這是1958年由Frank Rosenblatt所發展的感知機使用的輸出函數)作為轉換函數呢?

反向傳播及多層感知器